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ABSTRACT 

With the use of long-term numerical simulations, we study the evolution and orbital behavior 
of cometary nuclei in cold Kuiper belt-like debris disks under the gravitational influence of dwarf 
planets (DPs); we carry out these simulations with and without the presence of a Neptune-like 
giant planet. This exploratory study shows that in the absence of a giant planet, 10 DPs are 
enough to induce strong radial and vertical heating on the orbits of belt particles. On the other 
hand, the presence of a giant planet close to the debris disk, acts as a stability agent reducing 
the radial and vertical heating. With enough DPs, even in the presence of a Neptune-like giant 
planet some radial heating remains; this heating grows steadily, re-filling resonances otherwise 
empty of cometary nuclei. Specifically for the solar system, this secular process seems to be able 
to provide material that, through resonant chaotic diffusion, increase the rate of new comets 
spiraling into the inner planetary system, but only if more than the ~ 10 known DP sized objects 
exist in the trans-Neptunian region. 

Subject headings: planets and satellites: dynamical evolution and stability — Kuiper belt: general — 
methods: numerical 


1. Introduction 


One of the characteristics of evolved planetary 
systems is the prolonged presence of the rem¬ 
nants of stellar and planetary formation, rang¬ 
ing in size from dust grains to cometary nuclei 
to DPs. This material, located beyond the region 
where planets rapidly “clean-up” their vici nity, is 
known as a debris disk (for a review see IWvattI 
20081 : Kenvon et al.ll2008l and references therein). 


In our solar system the present day remnants in 
this region constitute the “Kuiper Belt” (KB). Al¬ 
though the lifetime of debris disks depends on di¬ 
verse factors, such as the stellar mass and neigh¬ 
boring environment, the majority of 100 Myr old 
stars have observational features consistent with 
the presence of debris disks and even a few 10 
Gyr old stars show evidence of having d ebris disks 
( Decin et al. 20031 : Greaves et al. 2005). 
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The first discovered extrasolar debris disk was 
the one of Vega, detected b y its infrared (IR) ex - 
cess with the IRAS satellite ( Aumann et al.l[l984h . 
The IR excess is believed to be produced by belts 
of dust particl es originating from a steady colli- 
sional cascade ( Muller et al. 201Clll : for the case of 
Vega, this belt is located at 100 AU from 
the star. The study of extrasolar debris disks is 
relevant in several aspects to the understanding of 
the planetary system formation process; moreover, 
debris disks have been employed to determine the 
presence of planets in extr asolar planetary systems 
( Zuckerman fc Son3l2004l l. 

On the other hand, DPs have an important role 
on the dynamics of primigenious planetary disks 
as the initiators of collisional cascades once they 
reach ^ 1000 to 3000 km size; they stir the orbits 
of residual icy planetesimals, increasing collisions; 
these collisions are responsible for both grounding 
some icy-planetesimals to dus t, as well as creating 
some super -Earth sized cores ( Kenvon &: BromlevI 
2004l20Rjl . Also, massive planets in evolved de¬ 
bris d isks are able to produce gaps an d dust out¬ 
flows ( Moro-Martin fc Malhotra 200, 'll) . 

In the specific case of the KB, recent stud¬ 
ies show that a number of its dynamical com¬ 
ponents ca n be explained with a migrating Nep¬ 
tune (e.g. Malhotral 1993: Levispn_^^Morbidelli 


2003 : Chiang et al.l l2007t Morbidelli et all 2008t 
Nesyorn-TT 2015ll . Indeed, all populations in the 
KB conserve evidences of their close interaction 
with the giant, except probably for the classical 
KB (CKB). The CKB has been defined as a bi- 
modal orbital distribution: the hot (incli nations 
i > 5°) and cold {i < 5°) components ( BrownI 
2 OOIII . However, some mixing bet ween both popu¬ 


lations seem to have taken place (iMorbidelli et al 
20081 : Volk fc Malhotral 2011 : Petit et ^ 2011 1. 


The most accepted scenario to explain the 
coexistence of both hot and cold p opulations 


coexistence or botn not and cold population; 
(iBatvgin et al. 2011 : Wolff et al. 2012t Nesvornv 


2015f) involves the action of a migrating Neptune, 
going outwards launching lots of planetesimals to 
form the hot population; the cold disk bodies, 
starting beyond 40 AU, simply kept their pri¬ 
mordial orbits mostly unaffected by Neptune that 
stopped migrating at some point in the evolution 
of the early solar system when the disk material 
grew scarce ( Gomes et al.ll2003l . 

Regarding the largest bodies of the power spec¬ 


tra on debris disks, the only examples we know 
are the KB objects (KBOs) with radii between 
400 and 1200 km, a fe w of which have on ly re¬ 


cently been discovered (iBrown et al.ll2005ll . Ex¬ 
trapolation of the size distribution of smaller 
KBOs has sometimes been used to attempt to 
estimate the numbers of such larger objects (i.e. 


Bernstein et al.l 1200411 . but estimations are still 


inconclusive. 

Regardless of their number, it is usually be¬ 
lieved DPs to have only a small influe nce on the 
evolut ion of debris disks in general. iFernand^ 
( 1980l l presents a first approximation where he at¬ 
taches importance to massive objects, of up to 
1.7 X lO-^Me, in a very massive KB disk (about 9 
Mq), finding that, in the presence of thousands of 
Ceres-like objects, direct encounters of cometary 
nuclei with larger bodies could lead to scatter of 
comets, sending them to the inner planetary re¬ 
gion, in this way possibly maintaining a steady 
influx of short-period comets. Current estimates 
of the mass and composition of the KB rule out 
this possibility as the main driver to produce the 
observed population of short-period comets. The 
infall inrate of comets on planetary systems might 
be of great importance in terms of habitability for 
example: it is believed that a large fraction of the 
water in t he primeval Earth ca me from comets and 
asteroids ( Altwegg et al. 201511 : also, at some later 
point it becomes necessary, for long-term evolu¬ 
tion of life, to have a reduced cometary infall rate. 
However at present, other than the KB, we are not 
able to observe such details on other debris disks. 


In this work we produce an exploratory study, 
that helps to better understand the dynamical ef¬ 
fects of DPs on cold Kuiper belt-like debris disks 
(KBLDD) with and without the influence of a 
Neptune-like giant planet. The physical effects 
presented here are of a general nature, as such, 
we expect them to be relevant in a wide variety of 
debris disks. In particular, we believe these results 
can be qualitatively applied to the KB (although 
we do not pretend to present a detailed study of 
the KB dynamics). A more quantitative study of 
the KB or of any other specific debris disk is be¬ 
yond the scope of this letter. 
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2. Simulations 


In this work we explore by means of long-term 
numerical simulations, the influence of random 
DPs on the dynamics of cold KBLDDs. The 
random DPs share physical characteristics with 
the ones observed in the solar system’s trans- 
Neptunian region, while the cold KBLDDs resem¬ 
ble the observed cold population of the solar sys¬ 
tem’s CKB. We constructed our initial conditions 
to resemble the cold CKB because it is the com¬ 
ponent least affected by Neptune, therefore the 
most stable. This is also the most intuitive start¬ 
ing point for a generic statistical study of debris 
disks. Among the differences with the solar sys¬ 
tem precise conditions are: the exact quantity of 
DPs, a zero inclination for our Neptune-like giant 
planet, and the random generated initial condi¬ 
tions of the belt particles. 


For our studies we employ the hybrid symplec- 
tic i ntegrator includ ed in the MERCURY pack¬ 
age ( Chamber^ 19991) . This integrator lets us fol¬ 
low the evolution of test particles in a potential 
generated by several major N-bodies plus a cen¬ 
tral star. It also permits to follow close encoun¬ 
ters between bodies with high accuracy by switch¬ 
ing from a symplectic to a Bulirsh-Stoer integra¬ 
tor; the switch between integrators takes place 
when particles get closer than a limit imposed 
in terms of the given major body’s Hill radius 
(Rh = (Mp/3Mq)1/3). 


All simulations are 1 Gyr long with an accuracy 
tolerance for the Bulirsh-Stoer integrator of 10“^°, 
a changeover distance between integrators of 3 Rh 
for any major body, and a time-step of 180 days for 
the symplectic integrator. The simulations were 
performed on 


2.1. Major Bodies 

The main central body in all simulations is a 1 
Mq star. 

We consider 4 different initial DP configura¬ 
tions: we use 10, 30, 50, and 100 randomly gener¬ 
ated cold DPs. The orbital parameters of all DPs 
lie within the following limits: semimajor axes, 
35 AU < a < 60 AU; eccentricities, 0.0 < e < 0.1; 
inclinations, 0.0° < i < 5.0°; arguments of peri- 


^ Atocatl is a supercomputer of the Institute de Astronomia 
at UNAM. 


center, 0° < w < 360°; longitudes of the ascending 
node, 0° < D < 360°; and mean anomalies, 0° < 
M < 360°. DP masses take random values in the 
range 3.3 x 1 O“®M 0 <m < 2.8 x IO^^Mq, where 
upper limit corresponds to Eris’s mass, while the 
lightest corresponds to the mass of 2002 AW 197 , 
this is, the biggest and one small but significative 
object in our trans-Neptunian region. 

All four DP configurations were run with and 
without the presence of a giant planet. The pa¬ 
rameters for this body were exactly the ones the 
real Neptune has but with zero inclination for the 
sake of simplicity, because the giant planet de¬ 
fines the angular momentum of the system (i.e. 
this represents the natural reference system of the 
problem); had we chosen different planes for the 
giant planet and the KBLDD an initial rearrange¬ 
ment of test particles would have occurred to come 
into balance with the giant planet’s plane. 


To better see the cumulative effect, we con¬ 
structed the sets of DPs in such a way that the 
larger DP sets include all the DPs of the previous 
set, i.e. the set of 10 DPs is a subset of the one 
of 30 DPs, etc. The total mass in DPs for 10, 30, 
50, and 100 objects is respectively: 0.011, 0.032, 
0.063, and 0.131 M 0 ; for co mparison, the CKB es- 
timated mass is ^ O.OIM 0 (IBernstein et al.ll2004 : 


Eraser et al.ll2014 1. 


2.2. Test particles’ initial conditions: Ran¬ 
dom Cold KBLDD 


We generate a belt of 1000 test particles that 
resemble th e observed cold CKB po p ulation. Ac¬ 
cordin g to iKavelaars et al. (l2008h. Petit et al. 
( 201ll) . and Dawson fc Murrav-Clav ( 20121) . the 
current cold CKB have orbits with semi-major 
axes, 42.5AU < a < 44.5AU, but mainly around 
44 AU, with inclinations, i < 4°, and eccentrici¬ 
ties, e < 0.05, for most objects of the population. 


We assign the values of the orbital parameters 
of the particles as follows: for a we use a random 
Gaussian distribution with mean and standard de¬ 
viation: (a) = 44.0 AU, tJa = 1-5 AU. For e and uj 
we generate a point distribution in an XY plane 
where each coordinate gets random Gaussian val¬ 
ues with mean zero and standard deviations given 
by cr(g^ gy) = 0.03; each point represents a vector, 
e = (ex,ev), whose magnitude, |ej = \/ex + Sy, 
is the e of the particle; also, we define the an- 
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gle between e and the X axis, 4)^, as w, therefore 
Lo = (f)e = TavT^icY!ex)\ in this manner the ini¬ 
tial e distribution has mean and standard devia¬ 
tion: (e) = 0.037, CTe = 0.019, while w is randomly 
distributed between 0° and 360°. We follow an 
analogous procedure to obtain the i and distri¬ 
butions; in this case we generate coordinates with 
random Gaussian points with mean zero and stan¬ 
dard deviations given by — 1-2°; the re¬ 

sulting i follows a distribution with (i) = 1.52°, 
ffi = 0.80°; while il is randomly distributed be¬ 
tween 0° and 360°. Finally, for M we use random 
values between 0° and 360°. 

3. Results and Discussion 

Fig. [T] shows the initial and final distributions 
of test particle eccentricities (left panel) and in¬ 
clinations (right panel) in the simulations without 
a Neptune-like giant planet; the black line rep¬ 
resents the initial conditions, while the different 
shades of blue represent the final distributions for 
10, 30, 50, and 100 DPs. Analogously, Fig.[^shows 
the same distributions when, along with the DPs, 
a Neptune-like planet is included at 30.09 A.U. 

From Fig. [T] we see that both e and i shift to¬ 
ward larger values as the number of DPs increases; 
this is to be expected as more DPs will produce 
a larger number of close encounters with test par¬ 
ticles, resulting in larger dispersions of e and i. 
A striking difference between e and i distributions 
can be noted: while for e there are more disturbed 
particles as the number of DPs increases; for i 
there seems to exist a saturation limit, where no 
particles can be heated beyond ^ 11°, not even 
with 100 DPs, while the mean of the distribution 
remains near ^ 5° with 30, 50, and 100 DPs. The 
latter is result of the initial distribution of DPs; 
as they are cold, with maximum initial inclina¬ 
tions of 5°, they do not seem to be able to push 
test particle’s i far beyond this limit. With 10 
DPs there is less dynamical heating and this limit 
is not reached, remaining around 4°. 

An interesting effect occurs when a Neptune- 
like planet is included in the simulations: as seen 
in Fig. [21 scattering induced by 10 and 30 DPs 
is severely damped for both e and i distributions. 
Again, with increasing DPs number, scattering of 
particles becomes stronger, leading to a shift of the 
distributions to higher values of e and i. For 50 


and 100 DPs, damping is slightly less important 
and, although fewer in number, some particles can 
rise to values of 0.20 and 11° for e and z, respec¬ 
tively (values similar to the ones reached without 
a giant planet). Again, the mean values of the 
final i distributions grow with DPs number, but 
always remain below 5°; even with 100 DPs, the 
mean is ~ 4°. This implies, contrary to intuition, 
that a giant planet can act as a stabilizing agent, 
by helping to vertically bound particles in its grav¬ 
itational potential (see Fig. 11). Mechanisms that 
could be responsible for this effect are: a) a su- 
pression on the number of close encounters of the 
cometary nuclei with DPs induced by the giant 
planet; from our studies we find an opposite be¬ 
havior, i.e. the presence of a giant planet increases 
the number of collisions due to the higher disk den¬ 
sity produced by its presence, b) Resonances with 
the giant; in this case, mean motion resonances 
(MMR) in the plane of the disk produced by the 
giant have a strong influence very high above the 
disk plane, flattening considerably the disk; this 
phenomenon has been recently demonstra ted to 
occur in galactic disks ( Moreno et al.l 2015 1. how¬ 
ever the lack of filamentary structure on Fig. |3l 
may suggest this effect is not important, c) Reso¬ 
nances induced by the DPs on the cometary nuclei; 
in this case the giant planet breakes the phases of 
the particle-DP interaction preventing the more 
efficient resonant heating, d) A gravitational non¬ 
resonant origin based only on the vertical force 
excerted by the giant; on average the giant acts 
like a 30 AU ring that pulls the cometary nuclei 
towards the plane of the disk producing the dis¬ 
tinctive triangle-like shape seen in Figure |3l 


With enough DPs, the effect of very close en¬ 
counters with DPs will be able to overcome the 
stabilizing influence of the giant planet; clearly, 
there must be a limit on how far this stabilizing 
influence can be exerted, but in the radii we ex¬ 
plore, we do not reach it. In the presence of the 
giant, there are more close encounters due to the 
higher density; this may lead to more dust pro¬ 
duction in the disk than without the presence of 
the giant. 

The left panels of Fig. |4| show the evolution 
throughout the simulations of (e) and (z), while 
the right panels show Oe and respectively. 
Broad lines show the evolution produced by DPs 
without a giant planet, while thin lines correspond 
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to simulations that include a Neptune-like body. 
The top-left panel of Figure |4] shows how, in all 
8 cases, (e) increases almost monotonically; nat¬ 
urally, as the number of DPs increases, their ef¬ 
fect on the final (e) increases. The top-right panel 
shows a similar behavior for ae (note the different 
scale between panels). These results strengthen 
what we have seen in the previous figures: the 
increasing presence of minor bodies increasingly 
perturbs the test particles, both with and without 
a giant planet. 

The growing radial heating allows test particles 
to encounter resonances, replenishing them with 
cometary nuclei. This effect is clear in spite of 
the small number of test particles we employ in 
our simulations. This becomes relevant not only 
because of the inherently fascinating behavior of 
particles trapped into resonances, but also because 
it is generally assumed that, in advanced stages of 
debris disks, there are no more known mechanisms 
able to restock the material on resonant regions. 

We also find that several of those particles are 
effectively trapped by resonances with the giant 
planet increasing dramatically their eccentricities. 
This mechanism might work as a plausible secular 
process able to sustain a rate of new comets spi¬ 
raling into the inner planetary system (this rate 
has not been fully explained for the KB). 

By comparing the thin to the broad lines in the 
two bottom panels of Fig. SI we can see the stabi¬ 
lizing effect of a Neptune-like planet: without the 
giant planet {i) quickly grows to reach the 5° limit 
found before, when a Neptune-like body is present 
evolution is smooth and rising but slower; with a 
giant planet, 100 DPs are required to produce a 
similar effect to what 10 DPs were able to achieve 
without the giant. Also, without a giant planet, 
30 DPs are enough to get close to some sort of 
saturation point, and there is very little difference 
between the final values for {i) for 50 and 100 DPs; 
the saturation value seems to be similar to the DPs 
inclination initial distribution. 

A similar trend is observed in the Ui evolution: 
the maximum dispersion reached is about 2.1° for 
30, 50, and 100 DPs without a giant planet, while 
with the giant this limit is about 1.6°. The effect 
produced by 10 bodies without the giant planet, 
clearly seen in both (i) and ui, almost disappears 
in the presence of the giant planet. In our so¬ 
lar system around 10 objects comparable in size 


to Pluto have been discovered, if this is the total 
number of this kind of bodies, their effect on our 
KB would be hardly noticeable; however, there is 
the possibility that the total number of DPs is sev¬ 
eral times larger. 


4. Conclusions 


With the use of long-term, N-body numerical 
simulations we have studied the dynamical effect 
of DPs on a cold debris disk with and without the 
presence of a giant planet. 

In the absence of the giant, DPs require only 1 
Gyr to induce substantial vertical heating on ini¬ 
tially cold test particles; this process increases the 
inclinations up to a saturation value of the order 
of the highest initial DP inclinations, in our sim¬ 
ulations, 5°. Likewise, radial heating (eccentricity 
dispersion) increases rapidly, although in this case, 
saturation is not reached. 


On the other hand, in the presence of a 
Neptune-like giant planet, the contribution of the 
DPs to the general heating diminishes severely. 
The 5° inclination limit obtained without the gi¬ 
ant planet is no longer reached, not even with 
100 DPs; in this case, the giant planet acts as 
a stability agent, concerning particle inclinations 
specifically, reducing the vertical heating. Regard¬ 
ing the radial heating, albeit a reduction is also 
observed, significant heating remains and grows 
steadily in time. The gravitational influence of 
the giant planet prevents the particles from dis¬ 
persing, keeping a higher density on the disk; this 
may have important consequences on the rate of 
collisions and on dust production. 


Another consequence of the heating produced 
by DPs is a slow but constant secular radial migra¬ 
tion of particles in the belt; several of those par¬ 
ticles are eventually trapped in the giant planet’s 
MMRs where, through chaotic diffusion, they 
could become part of other dynamical fa milies 
(e.g. Centaurs: iTiscareno fc Malhotra 2009ll . 


The continuous replenishing of resonant re¬ 
gions with new cometary nuclei leads several par¬ 
ticles through a dynamical evolution process that 
produces close encounters with the giant planet. 
Those bodies contribute to the influx rate of new 
short-period comets that may become important 
from the point of view of habitability, however 
observations of this mechanism are not yet avail- 
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able for planetary systems other than our Kuiper 
belt. In the case of the solar system this mech¬ 
anism may contribute to the short-period comet 
influx r ate, in better a ccordance with observa¬ 


tions ( Emehyanenko et al.ll2005l : IVolk fc Malhotra 


2008 . 2 OI 3 II : this is assuming the possibility of 
the existence of more than ten DPs in the trans- 
Neptunian region. Moreover, if the formation of 
several tens of DPs in the outer regions of our 
solar system took place prior to the migration of 
Neptune, a vertically pre-heated debris disk could 
have been already present when Neptune reached 
its current location; such process would produce 
a soft mixing between: the cold CKB population, 
the hot CKB population, and the resonant objects 
(those swept during Neptune’s migration). 
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Fig. 1.— Initial and final distributions of test par¬ 
ticles. Left panel shows the initial e distribution 
of KBLDD particles (black line) and the final dis¬ 
tributions, after 1 Gyr evolution, when 10 (darker 
blue line), 30 (middle blue line), 50 (lighter blue 
line), and 100 (gray line) random DPs are present 
in the simulation. Right panel is the same but for 
i. 



afAU) 
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Fig. 2.— Same as Fig.[T]but including a Neptune- 
like giant planet at 30.09 A.U. 


Fig. 3.— Height over the plane (z) vs. semimayor 
axis (a). The upper panel shows the cometary 
nuclei initial conditions (blue) and final (red) for 
our run with 50 DPs and no giant planet. Lower 
panel shows the same but for the run with 50 DPs 
plus the giant. 













































Fig. 4.— Evolution of test particle averages and standard deviations for eccentricity and inclination: top-left 
panel shows (e), top-right panel shows Ue, bottom-left shows (i), and bottom-right shows (7^. Broad lines 
stand for simulations without a giant planet, while thin lines for simulations that include a Neptune-like 
planet. Color code is the same as in Fig. [TJ 
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